{"cells":[{"cell_type":"markdown","metadata":{"id":"4511797C6EC64354AD05020541A54D1C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"# Lecture 14: Hierarchial Models  \n\n## Instructor: 胡传鹏（博士）[Dr. Hu Chuan-Peng]  \n\n### 南京师范大学心理学院[School of Psychology, Nanjing Normal University]"},{"cell_type":"markdown","metadata":{"id":"3ECD682A9279431B9934F253337C2450","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 回顾：广义线性模型GLM"},{"cell_type":"markdown","metadata":{"id":"FFF930555D0F438A889DE73495264C79","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"**当因变量为离散变量**  \n\n比如，因变量为二分变量 (答题正确率)，其中1代表回答正确，0代表回答错误。  \n\n正确率为离散变量，服从伯努利(Bernoulli)分布。  \n\n![Image Name](https://cdn.kesci.com/upload/image/rloa62fn8w.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"C0518C2B352F470BA7612606662B8611","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"对于因变量为离散变量的情况，我们需要使用广义线性模型(Generalized linear model，GLM)。  \n\n其特点为：  \n- 分布簇 (dist)不再局限于正态分布，而是允许其他不同的分布，比如 $y \\sim Bernoulli(p)$。  \n- 需要 **链接函数$g()$** 将 $\\alpha + \\beta * x$  映射到 $p$所在的范围。  \n\n| 一般线性模型 | 广义线性模型 |  \n|---|---|  \n| $y \\sim Normal(\\mu, \\sigma)$ | $y \\sim dist(p)$ |  \n| $\\mu = \\alpha + \\beta *x$ | $p = g(z)$|  \n|  | $z = \\alpha + \\beta *x$|  \n"},{"cell_type":"markdown","metadata":{"id":"7AFC54E93A50400989010177BCDC44EA","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"因变量：平均学习成绩（100制）  \n自变量：身高（cm）  \n\n$y_{score} \\sim N(\\mu, \\sigma)$  \n\n$\\mu = \\alpha + \\beta * height$  \n\n$\\alpha \\sim N(0, 200)$  \n\n$\\beta \\sim N(\\mu = 0, \\sigma = 10)$  \n\n$\\sigma \\sim HN(0, 10)$  \n\n数据结构：全国 --> 不同学校 --> 个体  \n\ncomplete pooling:  全国 --> 个体  \n\nNo pooling: 不同学校 --> 个体"},{"cell_type":"markdown","metadata":{"id":"4E8AD8AA26024BAFB274A51F928CBB7F","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"链接函数的具体转化过程，以逻辑(logit)回归为例：  \n1. 令 $z = \\alpha + \\beta *x$，$\\mu$的范围为 $(-\\infty, +\\infty)$。  \n2. $p = g(z)$，其中 $g()$ 为链接函数，输出结果 $p$ 的范围为 $(0,1)$。  \n3.  最后将 $p$ 输入到分布函数中 $y \\sim Bernoulli(p)$。  \n\n\n![Image Name](https://cdn.kesci.com/upload/image/rloa6zyf5a.png?imageView2/0/w/600/h/600)"},{"cell_type":"markdown","metadata":{"id":"22A6C1AF47B94886BD9D42C2DB086780","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"## Part 1: Hierarchical Models"},{"cell_type":"markdown","metadata":{"id":"99E5EDBCECC9419FAE018E860B8779E4","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"本节的目的在与：了解如何使用贝叶斯分层模型 (Bayesian Hierarchical model)对分层结构的数据进行分析。  \n\n重点在于：  \n- 了解分层模型的基本概念：包括变化截距与变化斜率。  \n- 通过 PyMC3 分别拟合变化截距模型与变化斜率模型(包含变化截距)。  \n"},{"cell_type":"markdown","metadata":{"id":"AC700BFAEA1B4CAB93CD63F56D01A31D","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 数据的层级结构  \n#### 我们想用工龄预测某个拥有4个学院的学校员工的工资  \n##### 收集到的数据如下：  \n\n![Image Name](https://cdn.kesci.com/upload/image/rmcutay7yg.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"1A86DFCDB18147958BDC4950AE206287","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"\n![Image Name](https://cdn.kesci.com/upload/image/rmcuwtjhi2.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"02DC28233BA145DC9014E6BFBC937449","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"如果需要使用工作年限预测工资水平，则最直觉的模型是线性模型：  \n$y \\sim Normal(\\mu, \\sigma)$  \n$\\mu  = \\alpha + \\beta *x_{years}$"},{"cell_type":"markdown","metadata":{"id":"536B387551CC41E1A127A75AAD1BB615","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"由于该数据具有层级结构，个体的数据嵌套在不同的院系之下，而不同院系嵌在学校之下。  \n\n考虑到院之间会有区别，而院系之内的个体之间可能同质性较高，因此，忽略院系可能会带来不准确的预测。  \n\n对于这种具有层级结构的数据，我们可能有不同的做法。  \n"},{"cell_type":"markdown","metadata":{"id":"9DBEB26C79544BB58A1E848302E9CF8C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 方式一：非池化参数(no pooling)"},{"cell_type":"markdown","metadata":{"id":"815D17D730B0408B97C26296B24C9F0C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"这种做法的假定是：各部门之间完全没有关联，不同部门之间彼此独立。  \n\n从统计上讲，是假定各组之间的参数没有关系，或者说是完全异质。  \n\n在这种情况下，该模型等价于每个部门各有一个广义线性模型。"},{"cell_type":"markdown","metadata":{"id":"93BE72E123514E8FBB2695CCE5ED42FD","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"即：  \n\n部门1：$y \\sim Normal(\\mu_1, \\sigma)$  \n$\\mu_1  = \\alpha_{1} + \\beta_{1} *x_{years}$  \n\n部门2：$y \\sim Normal(\\mu_1, \\sigma)$  \n$\\mu_1  = \\alpha_2 + \\beta_{2} *x_{years}$  \n\n...  \n\n部门j：$y \\sim Normal(\\mu_j, \\sigma)$  \n$\\mu  = \\alpha_j + \\beta_j *x_{years}$  \n"},{"cell_type":"markdown","metadata":{"id":"DA3316F9858542E0808403C807F237BE","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"\n![Image Name](https://cdn.kesci.com/upload/image/rmcwluwquv.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"0AE8DA0B666D492E8A2496B6C0C3AFE2","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 方式二：完全池化参数(complete pooling)"},{"cell_type":"markdown","metadata":{"id":"E8C95D30F21B4DECBC7AF968E5BE20A5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"这种做法假定各院系之间完全没有差异；  \n\n从数据分析角度而言，是认为各院系间共享同一个参数；即所有院系是完全同质的。  \n\n在这种情况下，该模型等价于每个院系共享一个广义线性模型。"},{"cell_type":"markdown","metadata":{"id":"950FAB211C304686A11C9BD0D91AEF5C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"即：  \n\n$y \\sim Normal(\\mu, \\sigma)$  \n$\\mu  = \\alpha + \\beta *x_{years}$  \n"},{"cell_type":"markdown","metadata":{"id":"EE23536D5D5D452DACF9F4A1DFE1A81A","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"\n![Image Name](https://cdn.kesci.com/upload/image/rmcx1f4dlx.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"154E6F119E3145B688FBB910BA472621","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 方式三：部分池化参数(partial pooling)"},{"cell_type":"markdown","metadata":{"id":"A4FD5B2316FF47F2A22D6958ECD2AD51","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"在这种方式中，我们假定各院系间有一定的关系，即他们是嵌套在同一个总体（学校）中。  \n\n部分池化参数就是各组的参数并非完全相同，但却是来自于某个特定分布。  \n\n分层模型（Hierarchcial model）即是对参数进行部分池化。"},{"cell_type":"markdown","metadata":{"id":"69DE6D069B41434FA94512E5BEAA3DAB","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"例如，对于正态分布模型来说，层级模型包括两部分：  \n\nGroup level (组水平):  \n\n$y_{i,j} \\sim Normal(\\mu_{j}, \\sigma_{j})$  \n\nPopulation level  (总体水平):  \n\n$\\mu_{j} \\sim Normal(\\mu, \\sigma)$"},{"cell_type":"markdown","metadata":{"id":"733F02B6E2D94F4F9685C2326F9288E0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"\n![Image Name](https://cdn.kesci.com/upload/image/rmcyhbglbb.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"93194318C2C440D59FA8D6CDF0EB5E65","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"在线性模型中，模型参数通常会分布截距与斜率，因此，理论上，层级线性模型中的截距与斜率均可以分层进行估计。  \n\n在上述例子中：变化截距与变化斜率会随着不同部门变化，表达式如下：  \n\n![Image Name](https://cdn.kesci.com/upload/image/rmdh2ogd7u.png?imageView2/0/w/640/h/640)"},{"cell_type":"markdown","metadata":{"id":"0504590214BE4A54A9B91514853ECF58","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"在实际中，分层模型参数是否进行分层估计本身是灵活的，例如，在上述例子中，可能会有三种情况：  \n\n- 不同部门员工工资存在不同。  \n\t- 这被称为变化截距，截距的变化代表不同部门间员工平均工资的差异。  \n\t- ![Image Name](https://cdn.kesci.com/upload/image/rmanj9sw62.png?imageView2/0/w/640/h/640)  \n- 不同学校中，员工工作年限对工资的预测作用会发生变化。  \n\t- 这被称为变化斜率，斜率在不同部门间的变化代表，工作年限对工资的预测作用会在不同部门将发生变化。  \n\t- ![Image Name](https://cdn.kesci.com/upload/image/rmank74ytm.png?imageView2/0/w/640/h/640)  \n- 变化截距 + 变化斜率  \n\t- ![Image Name](https://cdn.kesci.com/upload/image/rmankk5jgh.png?imageView2/0/w/640/h/640)  \n\n图片来源：http://mfviz.com/hierarchical-models/"},{"cell_type":"markdown","metadata":{"id":"6428814AED24416B9427751BF3A41F11","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"结合上节课的GLM和这节课的HLM，若因变量为离散变量，数据为嵌套结构时，我们可以讲两者进行结合，这就是GLMM."},{"cell_type":"markdown","metadata":{"id":"F4A0C7C173A2487081F68BDFF59DF3CC","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"\n![Image Name](https://cdn.kesci.com/upload/image/rmeb2c6lqc.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"C4A9E862CDFE4ECBA299BFF11880F7BF","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"## Part 2: 示例1——调查数据"},{"cell_type":"code","metadata":{"collapsed":false,"id":"D37209798C05422B8C44C83A84F7BE72","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#加载需要使用的库\n%matplotlib inline\nimport numpy as np \nfrom scipy import stats\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport arviz as az\nimport pymc3 as pm\n\nnp.random.seed(123)  # 变化数种子，确保随后生成的变化数相同","outputs":[{"output_type":"stream","name":"stderr","text":"WARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n"}],"execution_count":1},{"cell_type":"markdown","metadata":{"id":"B90C99E95CFC4C8D88188BEE0D45809E","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"**具体实例**  \n\n基于泰国初等教育数据库，我们想探究，学生是否接受学前教育对之后学生是否留级的影响？"},{"cell_type":"markdown","metadata":{"id":"ECBC1B50F14740C6A368CD16C1121E93","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"数据来源：泰国初等教育的全国调查 (Raudenbush & Bhumirat，1992)。  \n- 数据中的每一行表示一名学生。  \n- 因变量 `REPEAT` 为二分类变量，表示学生在初等教育期间是否留级，1 = 留级，0 = 非留级。  \n- 自变量为学生是否受过学前教育 `PPED` (0 = 否，1 = 是)。  \n\n数据来源：https://github.com/MultiLevelAnalysis/Datasets-third-edition-Multilevel-book/blob/master/chapter%206/Thaieduc/thaieduc.sav"},{"cell_type":"markdown","metadata":{"id":"9586141F89E949CCA7519FAF49BB0106","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### Workflow  \n\n在回顾广义线性模型 (GLM)与介绍分层线性模型 (HLM)的基础知识后，我们回到泰国初等教育的实例，并通过PyMC完成的贝叶斯分层广义线性模型建模的全过程 (full workflow)。  \n\n![Image Name](https://cdn.kesci.com/upload/image/rkvikqg9q6.png?imageView2/0/w/650/h/650)"},{"cell_type":"markdown","metadata":{"id":"2DDCDC39B6304A4E8A4A91BF0153D1D7","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (1) 提出研究问题  \n\n使用泰国初等教育数据探究的研究问题为：  \n- 学生是否接受学前教育对学生留级的影响？  \n- 在考虑不同学校的分层结构时，接受学前教育对学生是否留级的影响有哪些变化？  \n- 分层模型是否能提供额外的信息？"},{"cell_type":"markdown","metadata":{"id":"A50FD7804D524F64B3880426D5EE27B8","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (2) 数据收集"},{"cell_type":"markdown","metadata":{"id":"021AEE61F31040CFAB4031D7F1F16BA8","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"数据来源：泰国初等教育的全国调查 (Raudenbush & Bhumirat，1992)。  \n- 数据中的每一行表示一名学生。  \n- 因变量 `REPEAT` 为二分类变量，表示学生在初等教育期间是否留级，1 = 留级，0 = 非留级。  \n- 自变量为学生是否受过学前教育 `PPED` (0 = 否，1 = 是)。  \n- 其他变量包括学生性别 `SEX`、学生所在的学校 `SCHOOLID`和学校平均社会经济地位 (SES)分数 `MSESC`。  \n\n数据来源：https://github.com/MultiLevelAnalysis/Datasets-third-edition-Multilevel-book/blob/master/chapter%206/Thaieduc/thaieduc.sav"},{"cell_type":"code","metadata":{"collapsed":false,"id":"9B6390DCB0CD445FB2CC4E3CAE0C3BC4","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 加载数据\ndata = pd.read_csv(\"/home/mw/input/data1317/thaieduc.csv\")\ndata.PPED = data.PPED.map({'yes':1,'no':0})","outputs":[],"execution_count":2},{"cell_type":"code","metadata":{"collapsed":false,"id":"1CBFCC68AC9B47BF8C4C4E609C1604E5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"data.head()","outputs":[{"output_type":"execute_result","data":{"text/plain":"   SCHOOLID   SEX  PPED  REPEAT  MSESC\n0     10103  girl     0       0   0.88\n1     10103  girl     0       0   0.88\n2     10103  girl     1       0   0.88\n3     10103  girl     1       0   0.88\n4     10103  girl     1       0   0.88","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>SCHOOLID</th>\n      <th>SEX</th>\n      <th>PPED</th>\n      <th>REPEAT</th>\n      <th>MSESC</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>1</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>1</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>1</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":3}],"execution_count":3},{"cell_type":"markdown","metadata":{"id":"B41F251C7E6C4B179E708A6627D29ACF","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 描述统计与探索性分析"},{"cell_type":"code","metadata":{"collapsed":false,"id":"5C971162EE684D67AB3735D8D2F77D4B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# groupby 可以按照某条件对数据进行分组，之后再使用 mean函数可以实现分组求均值。\ndata.groupby('PPED').REPEAT.mean() ","outputs":[{"output_type":"execute_result","data":{"text/plain":"PPED\n0    0.277778\n1    0.191729\nName: REPEAT, dtype: float64"},"metadata":{},"execution_count":4}],"execution_count":4},{"cell_type":"markdown","metadata":{"id":"CA9B899E8F8B4A9081979BA9E75670CD","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"可以发现，没有接受过学前教育学生的留级率 (27.8%)高于接受过学前教育的学生 (19.2%)。  \n\n可视化如下："},{"cell_type":"code","metadata":{"collapsed":false,"id":"89018BCF2B4547C6A29D6D28C4A20421","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"data.groupby('PPED').REPEAT.mean().plot.bar()\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/89018BCF2B4547C6A29D6D28C4A20421/rmrx2qbp5.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":5},{"cell_type":"markdown","metadata":{"id":"B3DF575B0A254A9084E2988E1EE3C627","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"对比不同学校下，学前教育对留级率的影响。  \n\n其中，每一条线代表一个学校。  \n\n可见，不同学校间的学前教育对留级率的影响存在差异。但整体存在接受学前教育后留级率下降的趋势。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"654F08E84B2549A592536348AAA7F18A","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 计算不同学校下，是否接受学前教育的留级率\ndata_plot = data.groupby([\"SCHOOLID\",'PPED']).REPEAT.mean() \n# 根据留级率进行绘图\nfor schoolID,data_i in data_plot.groupby(\"SCHOOLID\"):\n    data_i.plot.line(x=\"PPED\")\nplt.xlabel(\"pped\")             # 横坐标为是否接受学前教育\nplt.xticks([0,1],[\"no\",\"yes\"]) # 0为否，1为是\nplt.ylabel(r\"REPEAT \\%\")       # 纵坐标为留级率","outputs":[{"output_type":"execute_result","data":{"text/plain":"Text(0, 0.5, 'REPEAT \\\\%')"},"metadata":{},"execution_count":6},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/654F08E84B2549A592536348AAA7F18A/rmrx2qc4ob.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":6},{"cell_type":"markdown","metadata":{"id":"DC0349DC2FF54C37B89B53D50FF28A07","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"思路：  \nComplete pooling  \n$y_i = \\alpha + \\beta * x_i$  \n"},{"cell_type":"markdown","metadata":{"id":"E915D5A87FBE4A258D6CE5E55C3255FA","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"No pooling:  \n\n学校1:  \n\n$y_{i, 1} = \\alpha_1 + \\beta_1 * x_{i, 1}$  \n\n...  \n\n学校500  \n"},{"cell_type":"markdown","metadata":{"id":"157C754BB77E4942911745E9D46892B3","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (3) 选择模型"},{"cell_type":"markdown","metadata":{"id":"4A4FC954571045C3BDBF55DAB79BE99A","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"由于因变量(是否留级)是二分离散变量，因此我们选择使用基于伯努利(Bernoulli)分布的广义线性模型(Generalized linear model，GLM)进行模型拟合。  \n\n这里我们考虑三种可能的模型：  \n- 模型1：仅包含变化截距的模型。表示拟合不同学校的平均留级率。  \n- 模型2：变化截距模型。在模型1的基础上，考虑学前教育在level1 (group)的总体效应，这种总体效应也成为固定效应。  \n- 模型3：变化斜率模型。在模型2的基础上，考虑学前教育在level2 (population)的不同学校上的不同效应，这种总体效应也成为变化效应。变化截距也属于变化效应。"},{"cell_type":"markdown","metadata":{"id":"C3C272DACB05457DB5A5675534FEF7BE","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 模型1：仅包含变化截距的模型"},{"cell_type":"markdown","metadata":{"id":"F848E1E316FB4240ACB6332A0096A7D2","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"notebookId":"638db1ba8cb7440653b762d4","trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"$$  \ny_{(i,j)} \\sim Bernoulli(p_{(i,j)})  \n$$  \n该模型只拟合不同学校的平均留级率。  \n- 其中$i$为不同学生，$j$为不同学校。  \n\n$$  \np_{(i,j)} = g(\\alpha_j) \\\\  \n\\alpha_j \\sim Normal(\\mu_{\\alpha},\\sigma_{\\alpha})  \n$$  \n- $\\alpha_{j}$ 为截距，代表平均留级率。该解决随着不同学校j进行变化，因此称为变化截距。  \n- g()为链接函数  \n- 其中，$\\mu_{\\alpha}, \\sigma_{\\alpha}$ 为 level 1 group 的效应，而 $\\alpha_{j}$ 为 level 2 population 的效应。  \n- 注意，此时的模型没有考虑自变量 **学前教育** 的影响。"},{"cell_type":"markdown","metadata":{"id":"C69743455E2D4BCB98A7EE6025D25D47","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"notebookId":"638db1ba8cb7440653b762d4","trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"**先验**  \n\n此时模型参数 ($\\alpha$)的先验与 $\\mu_{\\alpha}, \\sigma_{\\alpha}$有关。而 $\\mu_{\\alpha}, \\sigma_{\\alpha}$ 也服从于不同的分布。比如：  \n\n$$  \n\\mu_{\\alpha} \\sim Normal(0,100) \\\\  \n\\sigma_{\\alpha} \\sim HalfNormal(50)  \n$$  \n- 其中 $\\mu_{\\alpha}, \\sigma_{\\alpha}$ 称为超参 (hyperparameters)，可以理解为参数 $\\alpha$ 的参数。  \n- $Normal(0,100),HalfNormal(50)$ 称为超先验 (hyperprior)，即先验的先验。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"AA795AA7A048484C9E5DB01FE5AB2360","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 将数据分层变量\"学校(school)\"转换为因子(factor)类型\nschool_idxs, school = pd.factorize(data.SCHOOLID)\n# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\ncoords = {\n    \"school\": school,\n    \"obs_id\": np.arange(len(school_idxs)),\n}","outputs":[],"execution_count":7},{"cell_type":"code","metadata":{"collapsed":false,"id":"3AC314C61EEE4BA2AA51BA723820DF75","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model1:\n    # 定义超先验 Hyperpriors\n    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)  # 对应上述公式\n    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n\n    # 定义先验\n    alpha = pm.Normal('alpha', mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n    \n    # 定义数据分层变量\"学校(school)\"\n    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n    \n    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n    p = pm.Deterministic(\"p\", alpha[school_idx])\n    \n    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。\n    # 也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")","outputs":[],"execution_count":8},{"cell_type":"markdown","metadata":{"id":"241D70B10068482C97CF526BDADC3C6C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"我们通过 PyMC3自带的函数 `model_to_graphviz()` 将模型可视化。  \n\n需要注意的是：  \n- 我们假设 $\\alpha$ (先验)服从一个正态分布 其均值为 *mu_alpha*，误差为 *sigma_alpha*。  \n- 并且*mu_alpha* 和 *sigma_alpha* 各自服从不同的分布，如下图，这些分布为 **超先验(hyperpriors)**。  \n- 定义超先验的目的在于约束学校变量带来的差异，这被称为 \"shrinkage\"，是分层模型 partial pooling的效果。更多详情请参考 https://mc-stan.org/rstanarm/articles/pooling.html"},{"cell_type":"code","metadata":{"collapsed":false,"id":"AF14D9D2EC4849A883B512281049B95A","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model1)","outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7ff934739090>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/AF14D9D2EC4849A883B512281049B95A/rmrx2tv7y8.svg\">"},"metadata":{},"execution_count":9}],"execution_count":9},{"cell_type":"markdown","metadata":{"id":"ABD8FB3812394911B8E7E4C29DF081F8","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 模型2：变化截距模型  \n\n在模型1的基础上，考虑学前教育在level1的总体效应。  \n$$  \ny_{(i,j)} \\sim Bernoulli(p_{(i,j)}) \\\\  \np_{(i,j)} = g(\\alpha_j + \\beta*\\text{x}_{i, j}) \\\\  \n\\alpha_j \\sim Normal(\\mu_{\\alpha},\\sigma_{\\alpha})  \n$$  \n- $x$为自变量是否接受学前教育。  \n- 需要注意的是，学前教育的效应 $\\beta$ 不随学校进行变化。这意味着它在不同学生和不同学校间保持固定，因此也称为**固定效应**。  \n\n\n**先验**  \n\n$$  \n\\mu_{\\alpha} \\sim Normal(0,100) \\\\  \n\\sigma_{\\alpha} \\sim HalfNormal(50) \\\\  \n\\beta \\sim Normal(0,100) \\\\  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"5ADCD8F40AD14D69871CB131D91147A0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model2:\n    # 定义超先验 Hyperpriors\n    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)\n    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n\n    # 定义先验\n    alpha = pm.Normal('alpha',mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n    \n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n   \n    #####################################\n    # 与模型1的区别在于定义固定效应 beta, 以及线性公式p的变化\n    #####################################\n    beta = pm.Normal('beta',mu=0,sd=100)\n    p = pm.Deterministic(\"p\", alpha[school_idx] + beta * x)\n    \n    # Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")","outputs":[],"execution_count":10},{"cell_type":"code","metadata":{"collapsed":false,"id":"0644B67C7BDE462498716CA263F1BE64","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model2)","outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7ff936864990>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/0644B67C7BDE462498716CA263F1BE64/rmrx2twxip.svg\">"},"metadata":{},"execution_count":11}],"execution_count":11},{"cell_type":"markdown","metadata":{"id":"7D07B5ED4C434E3CB948D9D77FF57EE5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 模型3：变化截距+斜率模型  \n\n$$  \ny_{(i,j)} \\sim Bernoulli(p_{(i,j)}) \\\\  \n\np_{(i,j)} = g(\\alpha_j + \\beta_j*\\text{x}_{i, j}) \\\\  \n\n\\alpha_j \\sim Normal(\\mu_{\\alpha},\\sigma_{\\alpha}) \\\\  \n\n\\beta_j \\sim Normal(\\mu_{\\beta},\\sigma_{\\beta})  \n$$  \n在模型2的基础上，考虑学前教育在level2的不同学校上的不同效  \n- 此时，$x$的效应 $\\beta_{j}$ 随不同的学校 j 进行变化。  \n- 这种总体效应也成为变化效应。其实变化截距也属于变化效应。  \n\n\n**先验**  \n\n$$  \n\\mu_{\\alpha} \\sim Normal(0,100) \\\\  \n\\sigma_{\\alpha} \\sim HalfNormal(50) \\\\  \n\\mu_{\\beta} \\sim Normal(0,100) \\\\  \n\\sigma_{\\beta} \\sim HalfNormal(50) \\\\  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"3C874B02848A4463B26BB57BC4DACB5E","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model3:\n    #####################################\n    # 与模型2的区别在于，额外定义效应beta的超先验 Hyperpriors\n    #####################################\n    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)\n    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n    mu_beta = pm.Normal(\"mu_beta\", mu=0.0, sigma=100)\n    sigma_beta = pm.HalfNormal(\"sigma_beta\", 50)\n    \n    # 定义先验\n    alpha = pm.Normal('alpha', mu=mu_alpha, sd=sigma_alpha, dims=\"school\")\n    beta = pm.Normal('beta', mu=mu_beta, sd=sigma_beta, dims=\"school\")\n    \n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n    \n    #####################################\n    # 与模型2的区别在于，注意此时的beta随着学校school_idx进行变化\n    #####################################\n    p = pm.Deterministic(\"p\", alpha[school_idx] + beta[school_idx] * x)\n    \n    # Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")","outputs":[],"execution_count":12},{"cell_type":"markdown","metadata":{"id":"25FFA6E328A44D3CB78A5A2001AB43D1","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"注意观察，每个变化效应的先验 (alpha 和 beta)都存在超先验 (hyperpriors)。  \n\n同样，beta参数的超先验约束了不同学校间学前教育效应的影响，可避免出现极端值。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"34AFF8F58FCA4AD99FDB479DF35D029C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model3)","outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7ff936826250>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/34AFF8F58FCA4AD99FDB479DF35D029C/rmrx2ueqqb.svg\">"},"metadata":{},"execution_count":13}],"execution_count":13},{"cell_type":"markdown","metadata":{"id":"ABF9DF34B9A245F9988F5DA9ED3F82B0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"这里我们先考虑最后一个模型 (模型3)，在后面模型比较的部分，再考虑其他两个模型。"},{"cell_type":"markdown","metadata":{"id":"AF30EA899E1246DEA74B5726FABD3D14","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (4)选择先验"},{"cell_type":"code","metadata":{"collapsed":false,"id":"BB04D74376DF4F5E8E4CE6916C6E30FB","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model3:\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=1000)","outputs":[],"execution_count":14},{"cell_type":"markdown","metadata":{"id":"35D1E3CDB6B6465898AF62763F4CAEB5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"首先查看 p 的先验分布。  \n\n注意，这里没有用链接函数对 p 进行转换，因此其范围在 -200到200之间，而不是 0-1之间。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"E29A8AAA09294CC3928D3FE3F8F41251","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_density(\n    {'p':prior_checks['p']}\n    )\nplt.show()","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/base.py:174: UserWarning: More chains (1000) than draws (338). Passed array should have shape (chains, draws, *shape)\n  UserWarning,\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/E29A8AAA09294CC3928D3FE3F8F41251/rmrx2v98ye.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":15},{"cell_type":"markdown","metadata":{"id":"1E4F9BB1F6CA403082D3058BF9439A06","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"其次，查看模型截距和斜率的先验。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"96F8C253290B48CC8B969188D9367AC0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_density(\n    {'alpha':prior_checks['alpha'],\n    'beta':prior_checks['beta']}\n    )\nplt.show()","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/base.py:174: UserWarning: More chains (1000) than draws (19). Passed array should have shape (chains, draws, *shape)\n  UserWarning,\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 993.6x331.2 with 2 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/96F8C253290B48CC8B969188D9367AC0/rmrx2vkt7t.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":16},{"cell_type":"markdown","metadata":{"id":"0EE232B3A41A4F459FBB711AC0E98C22","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"最后，别忘了查看超先验所对应超参数(hyperparameters)的先验分布。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"3FDC6F8A3FBA44189CD0176A4194F470","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_density(\n    {\n        'mu_alpha':prior_checks['mu_alpha'],\n        'sigma_alpha':prior_checks['sigma_alpha'],\n        'mu_beta':prior_checks['mu_beta'],\n        'sigma_beta':prior_checks['sigma_beta']\n    }\n    )\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 1490.4x662.4 with 6 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/3FDC6F8A3FBA44189CD0176A4194F470/rmrx2vis8n.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":17},{"cell_type":"markdown","metadata":{"id":"A8A95151CB104CB28791362E775BBDED","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (5) 拟合数据"},{"cell_type":"markdown","metadata":{"id":"9159E395A6B741A191A8DBBCFC448D41","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"拟合数据需要注意，虽然我们主要探究的自变量为\"是否接受学前教育\"，但是数据分层结构变量\"学校\"也是非常重要的，不要忽视了。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"AD78B64D2A84434EA973F8D058187FB2","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"data.head()","outputs":[{"output_type":"execute_result","data":{"text/plain":"   SCHOOLID   SEX  PPED  REPEAT  MSESC\n0     10103  girl     0       0   0.88\n1     10103  girl     0       0   0.88\n2     10103  girl     1       0   0.88\n3     10103  girl     1       0   0.88\n4     10103  girl     1       0   0.88","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>SCHOOLID</th>\n      <th>SEX</th>\n      <th>PPED</th>\n      <th>REPEAT</th>\n      <th>MSESC</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>1</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>1</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>10103</td>\n      <td>girl</td>\n      <td>1</td>\n      <td>0</td>\n      <td>0.88</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":18}],"execution_count":18},{"cell_type":"code","metadata":{"collapsed":false,"id":"C0C6E12AE1C146B597796FC5539F7358","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"data.groupby('SCHOOLID').REPEAT.mean().plot.bar()\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/C0C6E12AE1C146B597796FC5539F7358/rmrx2wef8g.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":19},{"cell_type":"code","metadata":{"collapsed":false,"id":"A84FE8CF2F1F4B509ABF60D4E16F60D0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 将数据分层变量\"学校(school)\"转换为因子(factor)类型\nschool_idxs, school = pd.factorize(data.SCHOOLID)\n# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\ncoords = {\n    \"school\": school,\n    \"obs_id\": np.arange(len(school_idxs)),\n}\n\nwith pm.Model(coords=coords) as model3:\n    # 定义level2学校的 Hyperpriors\n    mu_alpha = pm.Normal(\"mu_alpha\", mu=0.0, sigma=100)\n    sigma_alpha = pm.HalfNormal(\"sigma_alpha\", 50)\n    mu_beta = pm.Normal(\"mu_beta\", mu=0.0, sigma=100)\n    sigma_beta = pm.HalfNormal(\"sigma_beta\", 50)\n    # 定义先验\n    alpha = pm.Normal('alpha',mu=mu_alpha,sd=sigma_alpha, dims=\"school\")\n    beta = pm.Normal('beta',mu=mu_beta,sd=sigma_beta, dims=\"school\")\n    # x为自变量，是之前已经载入的数据\n    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n    p = pm.Deterministic(\"p\", alpha[school_idx] + beta[school_idx] * x)\n    # Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")","outputs":[],"execution_count":20},{"cell_type":"markdown","metadata":{"id":"5B473A9EEE004985A6E18C77B434DD42","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (6)采样过程诊断  \n\n如果使用MCMC对后验进行近似，则需要首先对MCMC过程进行评估。  \n\n* 是否收敛；  \n* 是否接近真实的后验。  \n\n对采样过程的评估我们会采用目视检查或rhat这个指标"},{"cell_type":"code","metadata":{"collapsed":false,"id":"64F254C4F9B6496F9F939179CB09CA3F","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model3:\n    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n    trace3 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha, sigma_beta, mu_beta, sigma_alpha, mu_alpha]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 26 seconds.\n"},{"output_type":"stream","name":"stdout","text":"0, dim: obs_id, 338 =? 338\n"},{"output_type":"stream","name":"stderr","text":"There were 48 divergences after tuning. Increase `target_accept` or reparameterize.\nThere were 64 divergences after tuning. Increase `target_accept` or reparameterize.\nThe estimated number of effective samples is smaller than 200 for some parameters.\n"}],"execution_count":21},{"cell_type":"code","metadata":{"collapsed":false,"id":"01761BA8CD0D4B87A2FBE9A701D2606C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_trace(trace3, var_names=['alpha','beta'])\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 864x288 with 4 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/01761BA8CD0D4B87A2FBE9A701D2606C/rmrx3t4ikt.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":22},{"cell_type":"markdown","metadata":{"id":"E4772F96662B4F2785C30CCEC1E574FA","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"注意，由于参数会随着学校进行变化，因此每一所学校(19所学校)都对应一个参数，如下："},{"cell_type":"code","metadata":{"collapsed":false,"id":"DF1236A66F2640B0ACF878F260085748","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace3, var_names=['alpha','beta'], kind=\"diagnostics\")","outputs":[{"output_type":"execute_result","data":{"text/plain":"           mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nmu_alpha -1.043  0.329  -1.693   -0.462      0.017    0.012     393.0   \nmu_beta  -0.550  0.409  -1.316    0.191      0.021    0.015     400.0   \nalpha[0] -1.289  0.563  -2.374   -0.325      0.026    0.018     495.0   \nalpha[1] -0.968  0.453  -1.836   -0.104      0.018    0.013     648.0   \nalpha[2] -1.073  0.515  -2.003   -0.011      0.019    0.014     768.0   \n...         ...    ...     ...      ...        ...      ...       ...   \np[333]   -1.744  0.500  -2.685   -0.791      0.012    0.009    1966.0   \np[334]   -1.744  0.500  -2.685   -0.791      0.012    0.009    1966.0   \np[335]   -1.744  0.500  -2.685   -0.791      0.012    0.009    1966.0   \np[336]   -1.744  0.500  -2.685   -0.791      0.012    0.009    1966.0   \np[337]   -1.744  0.500  -2.685   -0.791      0.012    0.009    1966.0   \n\n          ess_tail  r_hat  \nmu_alpha     594.0   1.01  \nmu_beta      639.0   1.01  \nalpha[0]    1034.0   1.01  \nalpha[1]     865.0   1.00  \nalpha[2]    1077.0   1.01  \n...            ...    ...  \np[333]      1389.0   1.00  \np[334]      1389.0   1.00  \np[335]      1389.0   1.00  \np[336]      1389.0   1.00  \np[337]      1389.0   1.00  \n\n[380 rows x 9 columns]","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>mu_alpha</th>\n      <td>-1.043</td>\n      <td>0.329</td>\n      <td>-1.693</td>\n      <td>-0.462</td>\n      <td>0.017</td>\n      <td>0.012</td>\n      <td>393.0</td>\n      <td>594.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>mu_beta</th>\n      <td>-0.550</td>\n      <td>0.409</td>\n      <td>-1.316</td>\n      <td>0.191</td>\n      <td>0.021</td>\n      <td>0.015</td>\n      <td>400.0</td>\n      <td>639.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>alpha[0]</th>\n      <td>-1.289</td>\n      <td>0.563</td>\n      <td>-2.374</td>\n      <td>-0.325</td>\n      <td>0.026</td>\n      <td>0.018</td>\n      <td>495.0</td>\n      <td>1034.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>alpha[1]</th>\n      <td>-0.968</td>\n      <td>0.453</td>\n      <td>-1.836</td>\n      <td>-0.104</td>\n      <td>0.018</td>\n      <td>0.013</td>\n      <td>648.0</td>\n      <td>865.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>alpha[2]</th>\n      <td>-1.073</td>\n      <td>0.515</td>\n      <td>-2.003</td>\n      <td>-0.011</td>\n      <td>0.019</td>\n      <td>0.014</td>\n      <td>768.0</td>\n      <td>1077.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>...</th>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n    </tr>\n    <tr>\n      <th>p[333]</th>\n      <td>-1.744</td>\n      <td>0.500</td>\n      <td>-2.685</td>\n      <td>-0.791</td>\n      <td>0.012</td>\n      <td>0.009</td>\n      <td>1966.0</td>\n      <td>1389.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>p[334]</th>\n      <td>-1.744</td>\n      <td>0.500</td>\n      <td>-2.685</td>\n      <td>-0.791</td>\n      <td>0.012</td>\n      <td>0.009</td>\n      <td>1966.0</td>\n      <td>1389.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>p[335]</th>\n      <td>-1.744</td>\n      <td>0.500</td>\n      <td>-2.685</td>\n      <td>-0.791</td>\n      <td>0.012</td>\n      <td>0.009</td>\n      <td>1966.0</td>\n      <td>1389.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>p[336]</th>\n      <td>-1.744</td>\n      <td>0.500</td>\n      <td>-2.685</td>\n      <td>-0.791</td>\n      <td>0.012</td>\n      <td>0.009</td>\n      <td>1966.0</td>\n      <td>1389.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>p[337]</th>\n      <td>-1.744</td>\n      <td>0.500</td>\n      <td>-2.685</td>\n      <td>-0.791</td>\n      <td>0.012</td>\n      <td>0.009</td>\n      <td>1966.0</td>\n      <td>1389.0</td>\n      <td>1.00</td>\n    </tr>\n  </tbody>\n</table>\n<p>380 rows × 9 columns</p>\n</div>"},"metadata":{},"execution_count":62}],"execution_count":62},{"cell_type":"markdown","metadata":{"id":"8D99FC8BCDD745C985FB3F7B80B5D705","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (7)模型诊断  \n\n在MCMC有效的前提下，需要继续检验模型是否能够较好地拟合数据。  \n\n我们会使用后验预测分布通过我们得到的参数生成一批模拟数据，并将其与真实数据进行对比。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"3004A841D7EA497DA573F78BA9896D2C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 后验预测分布的计算仍在容器中进行\nwith model3:\n    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace3.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace3, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:03&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":24},{"cell_type":"code","metadata":{"collapsed":false,"id":"A8A2C5BB868747CEAF76D8552E484F72","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 绘制后验预测分布\naz.plot_ppc(trace3)\nplt.show()","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n  fig.canvas.print_figure(bytes_io, **kw)\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/A8A2C5BB868747CEAF76D8552E484F72/rmrx489nhj.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":25},{"cell_type":"markdown","metadata":{"id":"87A3D1505FE44B158A7D74322E39E751","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (8)模型比较"},{"cell_type":"markdown","metadata":{"id":"E58BEE47CE3E46348E4AA03E25119100","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"前面的诊断过程，我们只考虑了模型3 (变化截距+斜率模型)。在模型比较阶段，我们可以同时比较三个模型：  \n- 模型1：仅包含变化截距的模型，拟合不同学校的平均留级率。  \n- 模型2：变化截距模型。在模型1的基础上，考虑学前教育在group level的总体效应。  \n- 模型3：变化截距+斜率模型。在模型2的基础上，考虑学前教育在不同学校上 (population level)的不同效应。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"7857BCAAA27247C0AD49677064E5D6A5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 对模型进行mcmc采样\nwith model1:\n    trace1 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)\nwith model2:\n    trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)\n# 模型3前面完成了采样，这里可以不用再采样了\n# with model3:\n#     trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [alpha, sigma_alpha, mu_alpha]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 9 seconds.\nThere was 1 divergence after tuning. Increase `target_accept` or reparameterize.\nThere were 2 divergences after tuning. Increase `target_accept` or reparameterize.\nThe number of effective samples is smaller than 25% for some parameters.\n"},{"output_type":"stream","name":"stdout","text":"0, dim: obs_id, 338 =? 338\n"},{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [beta, alpha, sigma_alpha, mu_alpha]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 11 seconds.\nThere were 19 divergences after tuning. Increase `target_accept` or reparameterize.\nThere were 11 divergences after tuning. Increase `target_accept` or reparameterize.\nThe number of effective samples is smaller than 10% for some parameters.\n"},{"output_type":"stream","name":"stdout","text":"0, dim: obs_id, 338 =? 338\n"}],"execution_count":26},{"cell_type":"code","metadata":{"collapsed":false,"id":"0814C1599F0C4169A05925BA32290E9B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 将三个模型的采样结果进行比较\ncompare_dict = {\n    \"仅包含截距的模型\": trace1, \n    \"变化截距模型\": trace2, \n    \"变化斜率和截距模型\": trace3\n    }\n# 选择loo方法进行比较\ncomp = az.compare(compare_dict, ic='loo')\ncomp","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n"},{"output_type":"execute_result","data":{"text/plain":"           rank         loo      p_loo     d_loo   weight         se  \\\n仅包含截距的模型      0 -171.909766  10.250333  0.000000  0.51679  10.142879   \n变化斜率和截距模型     1 -171.964265  13.164235  0.054498  0.48321  10.384417   \n变化截距模型        2 -172.196791  10.979380  0.287025  0.00000  10.171711   \n\n                dse  warning loo_scale  \n仅包含截距的模型   0.000000    False       log  \n变化斜率和截距模型  1.608943    False       log  \n变化截距模型     1.174107    False       log  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>loo</th>\n      <th>p_loo</th>\n      <th>d_loo</th>\n      <th>weight</th>\n      <th>se</th>\n      <th>dse</th>\n      <th>warning</th>\n      <th>loo_scale</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>仅包含截距的模型</th>\n      <td>0</td>\n      <td>-171.909766</td>\n      <td>10.250333</td>\n      <td>0.000000</td>\n      <td>0.51679</td>\n      <td>10.142879</td>\n      <td>0.000000</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>变化斜率和截距模型</th>\n      <td>1</td>\n      <td>-171.964265</td>\n      <td>13.164235</td>\n      <td>0.054498</td>\n      <td>0.48321</td>\n      <td>10.384417</td>\n      <td>1.608943</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>变化截距模型</th>\n      <td>2</td>\n      <td>-172.196791</td>\n      <td>10.979380</td>\n      <td>0.287025</td>\n      <td>0.00000</td>\n      <td>10.171711</td>\n      <td>1.174107</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":27}],"execution_count":27},{"cell_type":"markdown","metadata":{"id":"2B852D8630FD4B629C11F1BAB3EE8FA0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"结果显示， m3(变化斜率和截距)的拟合度好于其他两个模型，这表示，在不同学校间学前教育对留级率的影响是不同的。"},{"cell_type":"markdown","metadata":{"id":"4083FD2431164C2DBF6E9DD8B30A75B4","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (9)统计推断"},{"cell_type":"markdown","metadata":{"id":"33CBACE929CF4133BF82FEFE1DD26B14","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"通过模型比较可以发现，学校所带来的变化效应的影响。如果体现这种影响呐？  \n\n通过 arviz提供的函数 `plot_forest()` 可以可视化在不同学校间学前教育所带来的影响。  \n- 图中左边的编号是 学校ID。  \n- 94% HDI 是学前教育效应 $\\beta$ 在不同学校间的可信区间。  \n- 通过 94% HDI判断效应是否可信发现，学前教育效应在一些学校间 (比如10418, 20204)存在差异，而在一些学校 (比如10109， 20309)不存在差异。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"E4862DAF62D44BF5B7F9404FF12636A3","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"sort_ = trace3.posterior[\"beta\"].mean(dim=[\"chain\", \"draw\"])\nsort_ = sort_.sortby(sort_, ascending=True)\n\naz.plot_forest(\n    trace3,\n    var_names=[\"beta\"],\n    combined=True,\n    coords={\"school\": sort_[\"school\"]}\n)\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x424.8 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/E4862DAF62D44BF5B7F9404FF12636A3/rmrx4z69bf.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":28},{"cell_type":"markdown","metadata":{"id":"AC0E91BCF5664E96B305B271FFE05F12","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"同样，我们可以通过 `az.summary` 函数查看各参数的具体数值。这里只展示了 beta参数的数值，大家可自行查看 alpha参数的相关统计值。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"81316A37EC3547E88EBF50C7E714BEC8","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":true,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace3, var_names = [\"beta\"])","outputs":[{"output_type":"execute_result","data":{"text/plain":"           mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nbeta[0]  -0.874  0.702  -2.297    0.355      0.030    0.021     523.0   \nbeta[1]  -0.166  0.617  -1.319    0.979      0.032    0.023     445.0   \nbeta[2]  -0.617  0.582  -1.749    0.428      0.021    0.015     749.0   \nbeta[3]   0.067  0.617  -1.037    1.250      0.033    0.023     360.0   \nbeta[4]  -1.168  0.850  -2.837    0.281      0.048    0.034     318.0   \nbeta[5]  -0.536  0.840  -2.190    1.011      0.026    0.018     968.0   \nbeta[6]  -0.842  0.661  -2.122    0.322      0.026    0.018     615.0   \nbeta[7]  -1.034  0.727  -2.539    0.262      0.039    0.028     341.0   \nbeta[8]  -0.419  0.682  -1.702    0.874      0.022    0.016     921.0   \nbeta[9]  -0.787  0.703  -2.153    0.546      0.029    0.020     578.0   \nbeta[10] -0.917  0.684  -2.279    0.283      0.030    0.021     501.0   \nbeta[11] -0.458  0.644  -1.620    0.842      0.024    0.017     724.0   \nbeta[12] -0.700  0.623  -1.847    0.484      0.022    0.015     793.0   \nbeta[13] -0.161  0.553  -1.153    0.917      0.022    0.016     609.0   \nbeta[14] -0.076  0.598  -1.189    1.026      0.025    0.018     562.0   \nbeta[15] -0.323  0.749  -1.842    1.084      0.028    0.020     748.0   \nbeta[16] -0.363  0.590  -1.430    0.813      0.025    0.018     554.0   \nbeta[17] -0.508  0.620  -1.682    0.683      0.021    0.015     904.0   \nbeta[18] -0.619  0.587  -1.679    0.509      0.024    0.017     608.0   \n\n          ess_tail  r_hat  \nbeta[0]      727.0   1.00  \nbeta[1]      454.0   1.01  \nbeta[2]     1376.0   1.00  \nbeta[3]      450.0   1.01  \nbeta[4]      506.0   1.00  \nbeta[5]     1455.0   1.00  \nbeta[6]      997.0   1.01  \nbeta[7]      560.0   1.00  \nbeta[8]     1926.0   1.00  \nbeta[9]      622.0   1.01  \nbeta[10]    1152.0   1.01  \nbeta[11]    1536.0   1.00  \nbeta[12]    1147.0   1.00  \nbeta[13]    1316.0   1.01  \nbeta[14]    1355.0   1.01  \nbeta[15]    1003.0   1.01  \nbeta[16]     497.0   1.01  \nbeta[17]    1477.0   1.00  \nbeta[18]     966.0   1.01  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>beta[0]</th>\n      <td>-0.874</td>\n      <td>0.702</td>\n      <td>-2.297</td>\n      <td>0.355</td>\n      <td>0.030</td>\n      <td>0.021</td>\n      <td>523.0</td>\n      <td>727.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[1]</th>\n      <td>-0.166</td>\n      <td>0.617</td>\n      <td>-1.319</td>\n      <td>0.979</td>\n      <td>0.032</td>\n      <td>0.023</td>\n      <td>445.0</td>\n      <td>454.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[2]</th>\n      <td>-0.617</td>\n      <td>0.582</td>\n      <td>-1.749</td>\n      <td>0.428</td>\n      <td>0.021</td>\n      <td>0.015</td>\n      <td>749.0</td>\n      <td>1376.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[3]</th>\n      <td>0.067</td>\n      <td>0.617</td>\n      <td>-1.037</td>\n      <td>1.250</td>\n      <td>0.033</td>\n      <td>0.023</td>\n      <td>360.0</td>\n      <td>450.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[4]</th>\n      <td>-1.168</td>\n      <td>0.850</td>\n      <td>-2.837</td>\n      <td>0.281</td>\n      <td>0.048</td>\n      <td>0.034</td>\n      <td>318.0</td>\n      <td>506.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[5]</th>\n      <td>-0.536</td>\n      <td>0.840</td>\n      <td>-2.190</td>\n      <td>1.011</td>\n      <td>0.026</td>\n      <td>0.018</td>\n      <td>968.0</td>\n      <td>1455.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[6]</th>\n      <td>-0.842</td>\n      <td>0.661</td>\n      <td>-2.122</td>\n      <td>0.322</td>\n      <td>0.026</td>\n      <td>0.018</td>\n      <td>615.0</td>\n      <td>997.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[7]</th>\n      <td>-1.034</td>\n      <td>0.727</td>\n      <td>-2.539</td>\n      <td>0.262</td>\n      <td>0.039</td>\n      <td>0.028</td>\n      <td>341.0</td>\n      <td>560.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[8]</th>\n      <td>-0.419</td>\n      <td>0.682</td>\n      <td>-1.702</td>\n      <td>0.874</td>\n      <td>0.022</td>\n      <td>0.016</td>\n      <td>921.0</td>\n      <td>1926.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[9]</th>\n      <td>-0.787</td>\n      <td>0.703</td>\n      <td>-2.153</td>\n      <td>0.546</td>\n      <td>0.029</td>\n      <td>0.020</td>\n      <td>578.0</td>\n      <td>622.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[10]</th>\n      <td>-0.917</td>\n      <td>0.684</td>\n      <td>-2.279</td>\n      <td>0.283</td>\n      <td>0.030</td>\n      <td>0.021</td>\n      <td>501.0</td>\n      <td>1152.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[11]</th>\n      <td>-0.458</td>\n      <td>0.644</td>\n      <td>-1.620</td>\n      <td>0.842</td>\n      <td>0.024</td>\n      <td>0.017</td>\n      <td>724.0</td>\n      <td>1536.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[12]</th>\n      <td>-0.700</td>\n      <td>0.623</td>\n      <td>-1.847</td>\n      <td>0.484</td>\n      <td>0.022</td>\n      <td>0.015</td>\n      <td>793.0</td>\n      <td>1147.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[13]</th>\n      <td>-0.161</td>\n      <td>0.553</td>\n      <td>-1.153</td>\n      <td>0.917</td>\n      <td>0.022</td>\n      <td>0.016</td>\n      <td>609.0</td>\n      <td>1316.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[14]</th>\n      <td>-0.076</td>\n      <td>0.598</td>\n      <td>-1.189</td>\n      <td>1.026</td>\n      <td>0.025</td>\n      <td>0.018</td>\n      <td>562.0</td>\n      <td>1355.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[15]</th>\n      <td>-0.323</td>\n      <td>0.749</td>\n      <td>-1.842</td>\n      <td>1.084</td>\n      <td>0.028</td>\n      <td>0.020</td>\n      <td>748.0</td>\n      <td>1003.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[16]</th>\n      <td>-0.363</td>\n      <td>0.590</td>\n      <td>-1.430</td>\n      <td>0.813</td>\n      <td>0.025</td>\n      <td>0.018</td>\n      <td>554.0</td>\n      <td>497.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>beta[17]</th>\n      <td>-0.508</td>\n      <td>0.620</td>\n      <td>-1.682</td>\n      <td>0.683</td>\n      <td>0.021</td>\n      <td>0.015</td>\n      <td>904.0</td>\n      <td>1477.0</td>\n      <td>1.00</td>\n    </tr>\n    <tr>\n      <th>beta[18]</th>\n      <td>-0.619</td>\n      <td>0.587</td>\n      <td>-1.679</td>\n      <td>0.509</td>\n      <td>0.024</td>\n      <td>0.017</td>\n      <td>608.0</td>\n      <td>966.0</td>\n      <td>1.01</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":29}],"execution_count":29},{"cell_type":"markdown","metadata":{"id":"861CDC89EA19482493E9752FEDB8A6E9","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"前面的分析是针对 population level 参数，即变化效应。而我们一般想了解组层面 (group level)的效应是否存在，即固定效应。  \n\n组层面的效应对应的参数是 mu_alpha, mu_beta, sigma_alpha, sigma_beta。我们可以通过 `az.summary` 与 `az.plot_posterior` 来查看这些参数。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"6FA38BF78962462DB9C1EC8F07AFDBFE","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_posterior(trace3, var_names=['mu_alpha','mu_beta'])\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 993.6x331.2 with 2 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/6FA38BF78962462DB9C1EC8F07AFDBFE/rmrx50qwh6.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":30},{"cell_type":"code","metadata":{"collapsed":false,"id":"DF466A3B2FE84F0E96EF925FD05D097C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":true,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace3, var_names=['mu_alpha','mu_beta'])","outputs":[{"output_type":"execute_result","data":{"text/plain":"           mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nmu_alpha -1.043  0.329  -1.693   -0.462      0.017    0.012     393.0   \nmu_beta  -0.550  0.409  -1.316    0.191      0.021    0.015     400.0   \n\n          ess_tail  r_hat  \nmu_alpha     594.0   1.01  \nmu_beta      639.0   1.01  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>mu_alpha</th>\n      <td>-1.043</td>\n      <td>0.329</td>\n      <td>-1.693</td>\n      <td>-0.462</td>\n      <td>0.017</td>\n      <td>0.012</td>\n      <td>393.0</td>\n      <td>594.0</td>\n      <td>1.01</td>\n    </tr>\n    <tr>\n      <th>mu_beta</th>\n      <td>-0.550</td>\n      <td>0.409</td>\n      <td>-1.316</td>\n      <td>0.191</td>\n      <td>0.021</td>\n      <td>0.015</td>\n      <td>400.0</td>\n      <td>639.0</td>\n      <td>1.01</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":31}],"execution_count":31},{"cell_type":"markdown","metadata":{"id":"18E440E164E3474884E4A6A55E666821","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"whoops！ 别忘了广义线性模型中链接函数的存在，我们需要把上述的参数值转换到0-1的范围。  \n\n可以通过 inv.logit函数进行转换，即 $1 / (1 + exp(-\\theta))$。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"2D56F25157604171935459D803580BC3","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 模型预测值\np_PPED0 = 1 / (1 + np.exp(-trace3.posterior[\"mu_alpha\"].mean())).to_pandas()\np_PPED1 = 1 / (1 + np.exp(-(trace3.posterior[\"mu_alpha\"].mean()+trace3.posterior[\"mu_beta\"].mean()))).to_pandas()\nprint(\n    \"未接受学前教育学生的留级率(%) = \",p_PPED0, \n    \"\\n 接受学前教育学生的留级率(%) = \", p_PPED1,\n    \"\\n 学前教育对留级率的影响(%) = \", p_PPED0 - p_PPED1\n    )","outputs":[{"output_type":"stream","name":"stdout","text":"未接受学前教育学生的留级率(%) =  0.26056411967000837 \n 接受学前教育学生的留级率(%) =  0.1689203056292715 \n 学前教育对留级率的影响(%) =  0.09164381404073688\n"}],"execution_count":32},{"cell_type":"markdown","metadata":{"id":"ABC62D3695FF434B8F7427F038B9B99B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"最后，我们比较，模型3(变化截距+斜率)和模型2(变化截距模型)中学前教育组层面的效应是否相同。  \n\n可见，两个模型的估计非常类似。但是变化截距+斜率模型比变化截距模型更多地考虑了学前教育效应在不同学校间的差异。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"74AB724E31CC4755A5DA96BE1B5BE37B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"beta1 = 1 / (1 + np.exp(-(trace2.posterior[\"beta\"].mean()))).to_pandas()\nbeta2 = 1 / (1 + np.exp(-(trace3.posterior[\"mu_beta\"].mean()))).to_pandas()\nprint(\n    \"变化截距模型中学前教育的效应(%) = \", beta1, \n    \"\\n 变化截距+斜率模型中学前教育的效应(%) = \", beta2\n    )","outputs":[{"output_type":"stream","name":"stdout","text":"变化截距模型中学前教育的效应(%) =  0.38296613205275537 \n 变化截距+斜率模型中学前教育的效应(%) =  0.3658040768754983\n"}],"execution_count":33},{"cell_type":"markdown","metadata":{"id":"4459946661714DC6BAA951B5D141F027","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 练习  \n\n前面我们定义了各种分层模型，但我们没有演示只包含变化斜率的模型，如下图。  \n\n接下来的练习大家可以尝试以下定义只包含变化斜率的贝叶斯广义线性分层模型。  \n\n关键在于：  \n- 只允许斜率 beta 随学习变化。而截距 alpha 不变化。  \n- 如何定义分层模型的先验以及相应的超先验。  \n\n![Image Name](https://cdn.kesci.com/upload/image/rmank74ytm.png?imageView2/0/w/640/h/640)"},{"cell_type":"markdown","metadata":{"id":"CA86FAC3B8A4490A973EF59A69BE8E7E","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"结合上节课学习的GLM和HLM，若我们的数据因变量为离散数据，且数据为嵌套结构时，我们可以将HLM和GLM结合。  \n\n"},{"cell_type":"code","metadata":{"collapsed":false,"id":"ECD911A9901E472FB3622F3426A05958","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#加载需要使用的库\n%matplotlib inline\nimport numpy as np \nfrom scipy import stats\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport arviz as az\nimport pymc3 as pm\n\nnp.random.seed(123)  # 变化数种子，确保随后生成的变化数相同\n\n# 加载数据\ndata = pd.read_csv(\"/home/mw/input/data1317/thaieduc.csv\")\ndata.PPED = data.PPED.map({'yes':1,'no':0})\n\n# 将数据分层变量\"学校(school)\"转换为因子(factor)类型\nschool_idxs, school = pd.factorize(data.SCHOOLID)\n# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\ncoords = {\n    \"school\": school,\n    \"obs_id\": np.arange(len(school_idxs)),\n}","outputs":[],"execution_count":34},{"cell_type":"code","metadata":{"collapsed":false,"id":"53366CB729664FD09016CF039F562FB7","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":true,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model_practice:\n    ############################################################################\n    # 要求：定义level2学校的 Hyperpriors\n    # 提示：通过 mu_beta = pm.Normal(???) 进行定义\n    ############################################################################\n    mu_beta = pm.Normal(\"mu_beta\", ...) \n    sigma_beta = pm.HalfNormal(\"sigma_beta\", ...)  # 注意这是 HalfNormal 分布\n    \n    ############################################################################\n    # 要求：定义先验\n    # 提示：通过 beta = pm.Normal(???) 进行定义。通过 dims=\"school\" 使得 beta随学校变化\n    ############################################################################\n    beta = pm.Normal('beta', ...)\n    alpha = pm.Normal('alpha', ...)\n\n    # 定义数据分层变量\"学校(school)\"\n    school_idx = pm.Data(\"school_idx\", school_idxs, dims=\"obs_id\")\n    x = pm.Data(\"x\", data.PPED, dims=\"obs_id\")\n    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n    p = pm.Deterministic(\"p\", alpha + beta[school_idx]*x)\n    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    y_obs = pm.Bernoulli(\"y_obs\", logit_p=p, observed=data.REPEAT, dims=\"obs_id\")","outputs":[{"output_type":"error","ename":"TypeError","evalue":"float() argument must be a string or a number, not 'ellipsis'","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/pymc3/theanof.py\u001b[0m in \u001b[0;36mfloatX\u001b[0;34m(X)\u001b[0m\n\u001b[1;32m     82\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 83\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloatX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     84\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mAttributeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mAttributeError\u001b[0m: 'ellipsis' object has no attribute 'astype'","\nDuring handling of the above exception, another exception occurred:\n","\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_1127/425640524.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;31m# 提示：通过 mu_beta = pm.Normal(???) 进行定义\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;31m############################################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m     \u001b[0mmu_beta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mNormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"mu_beta\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m...\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m     \u001b[0msigma_beta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mHalfNormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"sigma_beta\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m...\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# 注意这是 HalfNormal 分布\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/pymc3/distributions/distribution.py\u001b[0m in \u001b[0;36m__new__\u001b[0;34m(cls, name, *args, **kwargs)\u001b[0m\n\u001b[1;32m    119\u001b[0m             \u001b[0mdist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    120\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 121\u001b[0;31m             \u001b[0mdist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    122\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mVar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtotal_size\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdims\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    123\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/pymc3/distributions/distribution.py\u001b[0m in \u001b[0;36mdist\u001b[0;34m(cls, *args, **kwargs)\u001b[0m\n\u001b[1;32m    128\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mdist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    129\u001b[0m         \u001b[0mdist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__new__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 130\u001b[0;31m         \u001b[0mdist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    131\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mdist\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    132\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/pymc3/distributions/continuous.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, mu, sigma, tau, sd, **kwargs)\u001b[0m\n\u001b[1;32m    485\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtau\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_tensor_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtau\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    486\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 487\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmedian\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_tensor_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfloatX\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    488\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariance\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1.0\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtau\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    489\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/pymc3/theanof.py\u001b[0m in \u001b[0;36mfloatX\u001b[0;34m(X)\u001b[0m\n\u001b[1;32m     84\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mAttributeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     85\u001b[0m         \u001b[0;31m# Scalar passed\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 86\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloatX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     88\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/numpy/core/_asarray.py\u001b[0m in \u001b[0;36masarray\u001b[0;34m(a, dtype, order, like)\u001b[0m\n\u001b[1;32m    100\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_asarray_with_like\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morder\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlike\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlike\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 102\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    103\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    104\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mTypeError\u001b[0m: float() argument must be a string or a number, not 'ellipsis'"]}],"execution_count":35},{"cell_type":"code","metadata":{"collapsed":false,"id":"E7145A3ED8044289931808DBC0B7FC17","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model_practice)","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{"collapsed":false,"id":"AACBC3190B8A4009A669D16BF1524D4C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model_practice:\n    trace_practice = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{"collapsed":false,"id":"E35CB67133514F9D935B435D76C63074","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace_practice, var_names=['mu_beta'])","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{"collapsed":false,"id":"E9637BF787994FEDB08A834151E05E5B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_posterior(trace_practice, var_names=['mu_beta'])\nplt.show()","outputs":[],"execution_count":null},{"cell_type":"markdown","metadata":{"id":"978BD19C04034C4396857F564857E26B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"## Part 3: 示例2——信号检测论"},{"cell_type":"markdown","metadata":{"id":"1B596FB2B4F84165983EE49A81AB6B14","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### Workflow  \n\n![Image Name](https://cdn.kesci.com/upload/image/rkvikqg9q6.png?imageView2/0/w/650/h/650)"},{"cell_type":"markdown","metadata":{"id":"500A9B61534742C39C5B2B1C3440E539","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (1) 提出研究问题"},{"cell_type":"markdown","metadata":{"id":"F80EDBEBC5114DD7A7BBA67C2224FAE0","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"假如我们在进行一项再认实验，这个实验向被试呈现了一组图片，这些图片中有些是新的（被试之前没有见过的），而有些是旧的（被试之前已看过的）。  \n\n被试需要对呈现的每一张图片进行反应：“新（未见过）”或者“旧（见过）”。"},{"cell_type":"markdown","metadata":{"id":"4BF488D0FA144B67916052E877770237","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"信号检测论(Signal Detection Theory, SDT)通过分别测量被试的敏感性以及判断标准来建立每个被试再认过程的模型。  \n\n信号检测论模型的基本理念是：在每个试次中，当刺激呈现时，被试心里会产生一种“熟悉”（或者记忆强度）的信号。随后被试会依据这些信号，判断当前的刺激是“新”或“旧”。  \n![Image Name](https://th.bing.com/th/id/R.1565c09f95a7411070316c78242814c8?rik=rKBAJumpMrZ5Qg&riu=http%3a%2f%2fmedia.open.com.cn%2fmedia_file004%2f1003%2fdongshi%2fshiyanxlx%2fchapter6%2fimage%2f01.jpg&ehk=ui4XLQbTaXgMxP3BuNoBRAnww2cCrMLys2tN%2f1vm%2fvA%3d&risl=&pid=ImgRaw&r=0)  \n\n被试的反应有信号检测论的参数d'（辨别力指标/感觉敏感性指标）和 c（报告标准）决定。  \n"},{"cell_type":"markdown","metadata":{"id":"DB5942D815D5451DB841E48FCE6AD0D5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"那么我们该如何计算被试的参数d'（辨别力指标/感觉敏感性指标）和 c（报告标准）？"},{"cell_type":"markdown","metadata":{"id":"576E358FCEA3465DAAE7F87D252005AE","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (2) 数据收集"},{"cell_type":"markdown","metadata":{"id":"03813B5CEFF740B9ADA45BCFDB7ACD52","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"我们使用Skagerberg和Wright记忆一致性(memory conformity)研究中的对照组(Skagerberg & Wright，2008)"},{"cell_type":"code","metadata":{"collapsed":false,"id":"C91169BCEC35437C8963D1E9E0B2A72B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 加载数据\ndata = pd.read_csv(\"/home/mw/input/data1317/confcontr.csv\")\n# 查看数据\ndata.head()\n# 该数据中每一行都是一个试次，subno为被试编号，sayold为被试的新旧判断，isold为图片本身的新旧程度","outputs":[{"output_type":"execute_result","data":{"text/plain":"   Unnamed: 0  subno  sayold  isold\n0           1     53       1      0\n1           2     53       1      1\n2           3     53       1      1\n3           4     53       1      1\n4           5     53       1      0","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Unnamed: 0</th>\n      <th>subno</th>\n      <th>sayold</th>\n      <th>isold</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>1</td>\n      <td>53</td>\n      <td>1</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>2</td>\n      <td>53</td>\n      <td>1</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>3</td>\n      <td>53</td>\n      <td>1</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>4</td>\n      <td>53</td>\n      <td>1</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>5</td>\n      <td>53</td>\n      <td>1</td>\n      <td>0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":36}],"execution_count":36},{"cell_type":"markdown","metadata":{"id":"BCA48CDEA31C4CEDB6F9D47BF5305822","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 描述统计与探索性分析"},{"cell_type":"markdown","metadata":{"id":"D6852869A97044CCB4DBCA3B97D19D4F","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"首先，我们可以查看单个被试每个试次的反应（击中率(hit)、虚报率(false alarm)、正确拒斥率(correct rejection)、漏报率(miss)）"},{"cell_type":"code","metadata":{"collapsed":false,"id":"EE6411E72639497DBE4DF55FEA3493C5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"data = data.loc[data['subno']==53]\nhr = len(data.loc[(data['sayold']==1)&(data['isold']==1)])/len(data.loc[data['isold']==1])\nmr = len(data.loc[(data['sayold']==0)&(data['isold']==1)])/len(data.loc[data['isold']==1])\ncrr = len(data.loc[(data['sayold']==0)&(data['isold']==0)])/len(data.loc[data['isold']==0])\nfar = len(data.loc[(data['sayold']==1)&(data['isold']==0)])/len(data.loc[data['isold']==0])\nprint(f\"hr is {hr} \")\nprint(f\"mr is {mr} \")\nprint(f\"crr is {crr} \")\nprint(f\"far is {far} \")","outputs":[{"output_type":"stream","name":"stdout","text":"hr is 0.5319148936170213 \nmr is 0.46808510638297873 \ncrr is 0.6226415094339622 \nfar is 0.37735849056603776 \n"}],"execution_count":37},{"cell_type":"markdown","metadata":{"id":"7FB2E7EB21CD4C5CB0BE2A12DD2BF596","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (3) 选择模型"},{"cell_type":"markdown","metadata":{"id":"CFBB7FD677C040299E823B9EBF6FF82F","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"最常见的SDT模型，即假定被试对两种信号的熟悉程度的分布属于高斯分布（即正态分布），且这两个分布的方差是齐性的，但均值可能不同（例如，对旧刺激的熟悉程度更高）。"},{"cell_type":"markdown","metadata":{"id":"33A291E4A06948898101256E86AA2BA2","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"在方差齐性的情况下，对于单个被试，d'可由标准化后的击中率和虚报率(hit and false alarm rates)之差计算得来(Stanislaw & Todorov，1999)  \n\n$d'=\\phi^{-1}(HR)-\\phi^{-1}(FAR)$  \n\n$\\phi$是累积正态密度函数(cumulative normal density function)，可将z分数转换为概率，其反函数$\\phi^{-1}$ 将概率（例如击中率或虚报率）转换为z分数。标准化后的击中率和虚报率分别表示为 $z(HR)$ 和$z(FAR)$。  \n\n反应标准$c=-z(FAR)$（负虚报率的标准分数）表示(DeCarlo，1998)。  \n\n"},{"cell_type":"markdown","metadata":{"id":"ECEE574E925648CF9AECDC71262857B9","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"在除了直接对d'和c进行计算，我们可以使用广义线性模型来进行计算。  \n\n在二项反应的GLM中，我们假定反应属于Bernoulli分布.  \n\n即：$y_i \\sim Bernoulli(p_i)$"},{"cell_type":"markdown","metadata":{"id":"5DCBBD99ECBE400CA3E6C4A88E9CEAB9","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"被试反应即对图片新旧程度的判断是由击中率(HR)或虚报率(FAR)来决定的，故概率p表示击中率(HR)。  \n\n概率p本身的范围为0-1，因此，我们不能直接用线性模型直接对p进行建模。  \n\n相反，我们用一个连结函数将p转化为‘线性预测值’ $η$，并用线性回归模型对$η$建模。  \n\n即：$p_i = \\phi(\\eta_i)$  \n\n如前所述，$\\phi$是一个累积正态密度函数并能将z 值转化为概率，那我们可以用一个截距和斜率来表示$\\eta$。  \n"},{"cell_type":"markdown","metadata":{"id":"22E9DB53DA894CB99FC22EA3E08A9A5B","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"之所以用$\\phi$而不是logit函数，是因为击中率或虚报率本身是一个概率，它是由信号分布的z分数转化而来，而不是简单地将至于范围转化到[0,1]  \n\n即：$\\eta_i = \\beta_0 + \\beta_1 * isold_i$  \n\n由于$d'=\\phi^{-1}(HR)-\\phi^{-1}(FAR) = z(HR) - z(FAR)$，$c=-z(FAR)$  \n\n当isold为0，反应由虚报率决定，即$z(FAR) = -c$  \n\n当isold为1，反应由击中率决定，即$z(HR) = z(HR) - z(FAR) + z(FAR) = -c+d$  \n\n故在线性模型中，$\\beta_0$代表c(报告标准)的相反数，$\\beta_1$代表d(辨别力指标)  \n"},{"cell_type":"markdown","metadata":{"id":"629DCAA5EF084E42A102B7F3C6A6ECCA","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"我们首先展示一个单被试模型的例子，由于该模型只有一个被试，故是一个单层的模型  \n$$  \ny_{(i)} \\sim Bernoulli(p_{(i)}) \\\\  \np_{(i)} = \\phi(\\eta) \\\\  \n\\eta = -c + d *x_{isold}  \n$$  \n\n\n**先验**  \n\n$$  \n\\ c \\sim Normal(0,1) \\\\  \n\\ d \\sim Normal(0,1) \\\\  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"2942B445538B40979188350B2E587283","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model0:\n    # 定义priors\n    c = pm.Normal(\"c\", mu=0.0, sigma=1)  \n    d = pm.Normal(\"d\", mu=0.0, sigma=1)\n    # 自变量x是之前已经载入的数据\n    x = pm.Data(\"x\", data['isold'])\n    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n    eta = pm.Deterministic(\"eta\", -c + d*x)\n    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    sayold = pm.Bernoulli(\"syaold\", logit_p=eta, observed=data['sayold'])","outputs":[],"execution_count":38},{"cell_type":"code","metadata":{"collapsed":false,"id":"536F5E60042544BEA6681F085A570B3B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model0)","outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7ff9357d4f10>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/536F5E60042544BEA6681F085A570B3B/rmrx6rtdsk.svg\">"},"metadata":{},"execution_count":39}],"execution_count":39},{"cell_type":"markdown","metadata":{"id":"29430CCD839C411291EC5C13E4007BBF","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"单被试模型仅为简单说明如何用GLM建立信号检测论，接下来的过程为了避免混淆将不对该模型做进一步的分析"},{"cell_type":"markdown","metadata":{"id":"15A314BD99A04FF7AA07CEFC07107B12","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"接下来我们将展示在分层结构中，我们如何利用GLMM建立信号检测论的模型"},{"cell_type":"code","metadata":{"collapsed":false,"id":"8C1BC6BEA1F64654961914698CD72CA7","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"data = pd.read_csv(\"/home/mw/input/data1317/confcontr.csv\")\ndata.head()","outputs":[{"output_type":"execute_result","data":{"text/plain":"   Unnamed: 0  subno  sayold  isold\n0           1     53       1      0\n1           2     53       1      1\n2           3     53       1      1\n3           4     53       1      1\n4           5     53       1      0","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Unnamed: 0</th>\n      <th>subno</th>\n      <th>sayold</th>\n      <th>isold</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>1</td>\n      <td>53</td>\n      <td>1</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>2</td>\n      <td>53</td>\n      <td>1</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>3</td>\n      <td>53</td>\n      <td>1</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>4</td>\n      <td>53</td>\n      <td>1</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>5</td>\n      <td>53</td>\n      <td>1</td>\n      <td>0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":40}],"execution_count":40},{"cell_type":"markdown","metadata":{"id":"CC72098077C449188543A5F013539D0D","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 模型1：变化截距模型  \n\n我们假定各被试的c(报告标准)是不同的  \n\n**固定效应**  \n$$  \ny_{i,j} \\sim Bernoulli(p_{i,j}) \\\\  \np_{i,j} = \\phi(\\eta_{i,j}) \\\\  \n\\eta_{i,j} = -c_{i,j} + d_{i,j} *x_{isold,i,j}  \n$$  \n\n\n**先验**  \n\n$$  \n\\ c_{i,j} \\sim Normal(\\mu_{c,j},\\sigma_{c,j}) \\\\  \n\\ d_{i,j} \\sim Normal(0,1) \\\\  \n$$  \n**随机效应**  \n$$  \n\\ \\mu_{c,j} \\sim Normal(0,1) \\\\  \n\\ \\sigma_{c,j} \\sim HalfNormal(1) \\\\  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"4A2046FEF66343759E5F80DF274A8947","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 将数据分层变量\"被试(subject)\"转换为因子(factor)类型\nsubj_idxs, subject = pd.factorize(data.subno)\n# 定义学校与数据的映射：即标注哪名学生(行)属于哪一所学校\ncoords = {\n    \"subject\": subject,\n    \"obs_id\": np.arange(len(subj_idxs)),\n}","outputs":[],"execution_count":41},{"cell_type":"code","metadata":{"collapsed":false,"id":"6FC323DCFC25418B828CE7C36A31E158","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model1:\n    # 定义level2被试的 Hyperpriorspriors\n    mu_c = pm.Normal(\"mu_c\", mu=0.0, sigma=1)  \n    sigma_c = pm.HalfNormal(\"sigma_c\", sigma = 1)\n    # 定义priors\n    c = pm.Normal(\"c\", mu=mu_c, sigma=sigma_c,dims=\"subject\")  \n    d = pm.Normal(\"d\", mu=0.0, sigma=1)\n    # 定义数据分层变量\"被试(subject)\"\n    subj_idx = pm.Data(\"subject\", subj_idxs, dims=\"obs_id\")\n    # 自变量x是之前已经载入的数据\n    x = pm.Data(\"x\", data['isold'], dims=\"obs_id\")\n    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n    eta = pm.Deterministic(\"eta\", -c[subj_idx] + d*x)\n    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    sayold = pm.Bernoulli(\"syaold\", logit_p=eta, observed=data['sayold'],dims=\"obs_id\")","outputs":[],"execution_count":42},{"cell_type":"code","metadata":{"collapsed":false,"id":"2E61F69EDA584F4C9976C2943F9773E5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model1)","outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7ff7fd9af810>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/2E61F69EDA584F4C9976C2943F9773E5/rmrx6rvnkg.svg\">"},"metadata":{},"execution_count":43}],"execution_count":43},{"cell_type":"markdown","metadata":{"id":"1B6D5C54E1CE4D06B221DA325E35C7F5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"#### 模型2：变化斜率模型  \n\n我们假定各被试的d(辨别力指标)是不同的  \n\n**固定效应**  \n$$  \ny_{i,j} \\sim Bernoulli(p_{i,j}) \\\\  \np_{i,j} = \\phi(\\eta_{i,j}) \\\\  \n\\eta_{i,j} = -c_{i,j} + d_{i,j} *x_{isold,i,j}  \n$$  \n\n\n**先验**  \n\n$$  \n\\ c_{i,j} \\sim Normal(0,1) \\\\  \n\\ d_{i,j} \\sim Normal(\\mu_{d,j},\\sigma_{d,j}) \\\\  \n$$  \n**随机效应**  \n$$  \n\\ \\mu_{d,j} \\sim Normal(0,1) \\\\  \n\\ \\sigma_{d,j} \\sim HalfNormal(1) \\\\  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"70D753ADA6A7440A981BEE554C9268AD","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with pm.Model(coords=coords) as model2:\n    # 定义level2被试的 Hyperpriorspriors\n    mu_d = pm.Normal(\"mu_d\", mu=0.0, sigma=1)  \n    sigma_d = pm.HalfNormal(\"sigma_d\", sigma = 1)\n    # 定义priors\n    c = pm.Normal(\"c\", mu=0, sigma=1)  \n    d = pm.Normal(\"d\", mu=mu_d, sigma=sigma_d, dims=\"subject\")\n    # 定义数据分层变量\"被试(subject)\"\n    subj_idx = pm.Data(\"subject\", subj_idxs, dims=\"obs_id\")\n    # 自变量x是之前已经载入的数据\n    x = pm.Data(\"x\", data['isold'], dims=\"obs_id\")\n    # 定义线性模型：p是确定性变化变量，这个变量的值完全由右端值确定\n    eta = pm.Deterministic(\"eta\", -c + d[subj_idx]*x)\n    # 定义似然函数：Y的观测值，这是一个特殊的观测变化变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n    sayold = pm.Bernoulli(\"syaold\", logit_p=eta, observed=data['sayold'],dims=\"obs_id\")","outputs":[],"execution_count":44},{"cell_type":"code","metadata":{"collapsed":false,"id":"53FEFB17662E4319B98D9AB40098807B","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(model2)","outputs":[{"output_type":"execute_result","data":{"text/plain":"<graphviz.graphs.Digraph at 0x7ff7fd87c4d0>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/53FEFB17662E4319B98D9AB40098807B/rmrx6rwp5h.svg\">"},"metadata":{},"execution_count":45}],"execution_count":45},{"cell_type":"markdown","metadata":{"id":"FA9C6CA1A8DC4600881165ECB137889D","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (4)选择先验"},{"cell_type":"code","metadata":{"collapsed":false,"id":"61016CCDE10A426392CB4E07F6A9FE81","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model1:\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=1000)","outputs":[],"execution_count":46},{"cell_type":"code","metadata":{"collapsed":false,"id":"CCE7324A345448508A639B7181F2B557","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_density(\n    {'eta':prior_checks['eta']}\n    )\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/CCE7324A345448508A639B7181F2B557/rmrx6sk1pi.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":47},{"cell_type":"code","metadata":{"collapsed":false,"id":"973D6051ACE045339703FAE947F6E465","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model2:\n    # 先验预测检查\n    prior_checks = pm.sample_prior_predictive(samples=1000)","outputs":[],"execution_count":48},{"cell_type":"code","metadata":{"collapsed":false,"id":"15214C071D7944CD9EEF7833B5DAF992","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_density(\n    {'eta':prior_checks['eta']}\n    )\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/15214C071D7944CD9EEF7833B5DAF992/rmrx6tjxvi.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":49},{"cell_type":"markdown","metadata":{"id":"29A0EBF451A44567BD3A3A1052476C93","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (5) 拟合数据"},{"cell_type":"code","metadata":{"collapsed":false,"id":"6932A233D0614FAD912D64A44B32E809","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model1:\n    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n    trace1 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [d, c, sigma_c, mu_c]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 16 seconds.\n"},{"output_type":"stream","name":"stdout","text":"0, dim: obs_id, 3100 =? 3100\n"}],"execution_count":50},{"cell_type":"code","metadata":{"collapsed":false,"id":"A5FD3B2CF651489F9A74145672068269","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"with model2:\n    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n    trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)","outputs":[{"output_type":"stream","name":"stderr","text":"Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [d, c, sigma_d, mu_d]\n"},{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"stream","name":"stderr","text":"Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 15 seconds.\n"},{"output_type":"stream","name":"stdout","text":"0, dim: obs_id, 3100 =? 3100\n"}],"execution_count":51},{"cell_type":"markdown","metadata":{"id":"98534988945B4702847319B83FB6880D","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (6)采样过程诊断  \n\n如果使用MCMC对后验进行近似，则需要首先对MCMC过程进行评估。  \n\n* 是否收敛；  \n* 是否接近真实的后验。  \n\n对采样过程的评估我们会采用目视检查或rhat这个指标"},{"cell_type":"code","metadata":{"collapsed":false,"id":"9E188EE34B9243E5B344700F9E86137D","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace1, var_names = [\"mu_c\",\"sigma_c\",\"d\"])","outputs":[{"output_type":"execute_result","data":{"text/plain":"          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nmu_c     1.053  0.081   0.906    1.208      0.002    0.001    1638.0   \nsigma_c  0.332  0.066   0.214    0.456      0.001    0.001    1899.0   \nd        1.678  0.079   1.535    1.830      0.002    0.001    1537.0   \n\n         ess_tail  r_hat  \nmu_c       2564.0    1.0  \nsigma_c    2027.0    1.0  \nd          2739.0    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>mu_c</th>\n      <td>1.053</td>\n      <td>0.081</td>\n      <td>0.906</td>\n      <td>1.208</td>\n      <td>0.002</td>\n      <td>0.001</td>\n      <td>1638.0</td>\n      <td>2564.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>sigma_c</th>\n      <td>0.332</td>\n      <td>0.066</td>\n      <td>0.214</td>\n      <td>0.456</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>1899.0</td>\n      <td>2027.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>d</th>\n      <td>1.678</td>\n      <td>0.079</td>\n      <td>1.535</td>\n      <td>1.830</td>\n      <td>0.002</td>\n      <td>0.001</td>\n      <td>1537.0</td>\n      <td>2739.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":52}],"execution_count":52},{"cell_type":"code","metadata":{"collapsed":false,"id":"EE458F3390A241C290153269D0533410","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace2, var_names = [\"mu_d\",\"sigma_d\",\"c\"])","outputs":[{"output_type":"execute_result","data":{"text/plain":"          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\nmu_d     1.663  0.120   1.449    1.894      0.003    0.002    1772.0   \nsigma_d  0.493  0.101   0.320    0.695      0.002    0.002    1809.0   \nc        1.033  0.056   0.932    1.141      0.001    0.001    1703.0   \n\n         ess_tail  r_hat  \nmu_d       2216.0    1.0  \nsigma_d    1881.0    1.0  \nc          2430.0    1.0  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>mu_d</th>\n      <td>1.663</td>\n      <td>0.120</td>\n      <td>1.449</td>\n      <td>1.894</td>\n      <td>0.003</td>\n      <td>0.002</td>\n      <td>1772.0</td>\n      <td>2216.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>sigma_d</th>\n      <td>0.493</td>\n      <td>0.101</td>\n      <td>0.320</td>\n      <td>0.695</td>\n      <td>0.002</td>\n      <td>0.002</td>\n      <td>1809.0</td>\n      <td>1881.0</td>\n      <td>1.0</td>\n    </tr>\n    <tr>\n      <th>c</th>\n      <td>1.033</td>\n      <td>0.056</td>\n      <td>0.932</td>\n      <td>1.141</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>1703.0</td>\n      <td>2430.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":53}],"execution_count":53},{"cell_type":"code","metadata":{"collapsed":false,"id":"1C51C0570C844EBDA3721A5831D951B1","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_trace(trace1, var_names = [\"mu_c\",\"sigma_c\",\"d\"])\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 864x432 with 6 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/1C51C0570C844EBDA3721A5831D951B1/rmrx7ww58h.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":54},{"cell_type":"code","metadata":{"collapsed":false,"id":"77B5DDE264304E2389C31DF074AD82D9","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_trace(trace2, var_names = [\"mu_d\",\"sigma_d\",\"c\"])\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 864x432 with 6 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/77B5DDE264304E2389C31DF074AD82D9/rmrx7xieqa.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":55},{"cell_type":"markdown","metadata":{"id":"D54E9485BE2E42BC980036BADAE7814E","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (7)模型诊断  \n\n在MCMC有效的前提下，需要继续检验模型是否能够较好地拟合数据。  \n\n我们会使用后验预测分布通过我们得到的参数生成一批模拟数据，并将其与真实数据进行对比。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"FE593957A6794B8987A32D895A07C579","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 后验预测分布的计算仍在容器中进行\nwith model1:\n    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace1.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace1, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:04&lt;00:00]\n    </div>\n    "},"metadata":{}},{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"}],"execution_count":56},{"cell_type":"code","metadata":{"collapsed":false,"id":"205EC5BBCDBD41AA973F6A6C1B3C5F43","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 后验预测分布的计算仍在容器中进行\nwith model2:\n    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n    ppc_y = pm.sample_posterior_predictive(trace2.posterior) \n#将ppc_y转化为InferenceData对象合并到trace中\naz.concat(trace2, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)","outputs":[{"output_type":"display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n<style>\n    /* Turns off some styling */\n    progress {\n        /* gets rid of default border in Firefox and Opera. */\n        border: none;\n        /* Needs to be in here for Safari polyfill so background images work as expected. */\n        background-size: auto;\n    }\n    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n    }\n    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n        background: #F44336;\n    }\n</style>\n"},"metadata":{}},{"output_type":"update_display_data","data":{"text/plain":"<IPython.core.display.HTML object>","text/html":"\n    <div>\n      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      100.00% [4000/4000 00:03&lt;00:00]\n    </div>\n    "},"metadata":{}}],"execution_count":57},{"cell_type":"code","metadata":{"collapsed":false,"id":"5113C1D98E7A4141A8301166A11BB583","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 绘制后验预测分布\naz.plot_ppc(trace1)\nplt.show()","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n  fig.canvas.print_figure(bytes_io, **kw)\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/5113C1D98E7A4141A8301166A11BB583/rmrx8iq3ww.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":58},{"cell_type":"code","metadata":{"collapsed":false,"id":"C2D5C95EB601416E982EA925073D7EB1","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 绘制后验预测分布\naz.plot_ppc(trace2)\nplt.show()","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n  fig.canvas.print_figure(bytes_io, **kw)\n"},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/C2D5C95EB601416E982EA925073D7EB1/rmrx8rbrdw.png\">"},"metadata":{"needs_background":"light"}}],"execution_count":59},{"cell_type":"markdown","metadata":{"id":"13B8B5C8070D4372805D702737B1E7C5","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### (8)模型比较"},{"cell_type":"markdown","metadata":{"id":"60E4F949FBFD4ABDBEB4519CA79F4D9C","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"在模型比较阶段，我们可以同时比较2个模型：  \n- 模型1：变化截距模型。  \n- 模型2：变化截距+斜率模型。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"85FDABC4422341CF9A61C7B82DCA2C77","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 将2个模型的采样结果进行比较\ncompare_dict = {\n    \"变化截距模型\": trace1, \n    \"变化截距+斜率模型\": trace2\n    }\n# 选择loo方法进行比较\ncomp = az.compare(compare_dict, ic='loo')\ncomp","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n"},{"output_type":"execute_result","data":{"text/plain":"           rank          loo      p_loo     d_loo    weight         se  \\\n变化截距+斜率模型     0 -1869.178154  23.466945  0.000000  0.525785  22.262004   \n变化截距模型        1 -1869.965300  22.893202  0.787145  0.474215  22.222396   \n\n                dse  warning loo_scale  \n变化截距+斜率模型  0.000000    False       log  \n变化截距模型     6.057976    False       log  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>loo</th>\n      <th>p_loo</th>\n      <th>d_loo</th>\n      <th>weight</th>\n      <th>se</th>\n      <th>dse</th>\n      <th>warning</th>\n      <th>loo_scale</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>变化截距+斜率模型</th>\n      <td>0</td>\n      <td>-1869.178154</td>\n      <td>23.466945</td>\n      <td>0.000000</td>\n      <td>0.525785</td>\n      <td>22.262004</td>\n      <td>0.000000</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>变化截距模型</th>\n      <td>1</td>\n      <td>-1869.965300</td>\n      <td>22.893202</td>\n      <td>0.787145</td>\n      <td>0.474215</td>\n      <td>22.222396</td>\n      <td>6.057976</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":60}],"execution_count":60},{"cell_type":"code","metadata":{"collapsed":false,"id":"AF0A8FE561504E6DAC090F02F42FD2AB","jupyter":{},"notebookId":"638db1ba8cb7440653b762d4","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 将2个模型的采样结果进行比较\ncompare_dict = {\n    \"变化截距模型\": trace1, \n    \"变化截距+斜率模型\": trace2\n    }\n# 选择loo方法进行比较\ncomp = az.compare(compare_dict, ic='waic')\ncomp","outputs":[{"output_type":"stream","name":"stderr","text":"/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n  \"The default method used to estimate the weights for each model,\"\n"},{"output_type":"execute_result","data":{"text/plain":"           rank         waic     p_waic    d_waic    weight         se  \\\n变化截距+斜率模型     0 -1869.167259  23.456050  0.000000  0.526047  22.261703   \n变化截距模型        1 -1869.964040  22.891943  0.796781  0.473953  22.222379   \n\n                dse  warning waic_scale  \n变化截距+斜率模型  0.000000    False        log  \n变化截距模型     6.057708    False        log  ","text/html":"<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>rank</th>\n      <th>waic</th>\n      <th>p_waic</th>\n      <th>d_waic</th>\n      <th>weight</th>\n      <th>se</th>\n      <th>dse</th>\n      <th>warning</th>\n      <th>waic_scale</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>变化截距+斜率模型</th>\n      <td>0</td>\n      <td>-1869.167259</td>\n      <td>23.456050</td>\n      <td>0.000000</td>\n      <td>0.526047</td>\n      <td>22.261703</td>\n      <td>0.000000</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n    <tr>\n      <th>变化截距模型</th>\n      <td>1</td>\n      <td>-1869.964040</td>\n      <td>22.891943</td>\n      <td>0.796781</td>\n      <td>0.473953</td>\n      <td>22.222379</td>\n      <td>6.057708</td>\n      <td>False</td>\n      <td>log</td>\n    </tr>\n  </tbody>\n</table>\n</div>"},"metadata":{},"execution_count":61}],"execution_count":61},{"cell_type":"markdown","metadata":{"id":"8CA7D90095DB4FAA81B40C8A7AD5520D","notebookId":"638db1ba8cb7440653b762d4","jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"trusted":true,"runtime":{"status":"default","execution_status":null}},"source":"### 总结  \n\n- 本节课学习了分层线性模型的基本概念：包括变化截距与变化斜率、固定效应与变化效应、group level与population level。  \n- 通过 PyMc 分别拟合变化截距模型与变化斜率模型(包含变化截距)。"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python","nbconvert_exporter":"python","file_extension":".py","version":"3.5.2","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":0}